#!/usr/bin/python3

from pylab import *
import scipy.integrate


def conventional():
    """
    solve Cauchy's problem

      y'(t) = f(t, y)
      y_0 = y(t_0)

    """
    f = lambda t, y: cos(t)
    y_0 = 0.0
    t_max = 10.0
    dt = 0.05
    T = []
    Y = []

    integr = scipy.integrate.ode(f)
    integr.set_integrator('vode').set_initial_value(y_0, t=0.0)
    T.append(0.0)
    Y.append(y_0)

    while integr.successful() and integr.t < t_max:
        integr.integrate(integr.t + dt)
        T.append(integr.t)
        Y.append(integr.y)

    plot(T, Y, marker='.', label='conventional')

def recurrence_relation(h=0.01):
    '''
    unconventional
    '''
    f = lambda t, y: cos(t)
    T = []
    Y = []
    y_0 = 0.0
    t_0 = 0.0
    t_max = 10.0

    t = t_0
    y = y_0
    while t < t_max:
        T.append(t)
        Y.append(y)
        y = y + h * f(t, y)
        t = t + h

    plot(T, Y, marker=',', label='recurrence relation')

if __name__ == "__main__":
    D = arange(0, 10, 0.05)
    plot(D, sin(D), color='#CC00CC', linewidth=0.5, label='sin')
    #conventional()
    recurrence_relation(0.05)
    recurrence_relation(0.0001)
    xlim(min(D - 0.1), max(D + 0.1))
    ylim(min(sin(D)) * 2, max(sin(D)) * 2)
    legend()
    grid()
    show()
